###Function to estimate a panel model with a geospatial and temporal lag.###
###By Trey Hood and Jamie Monogan###

##Directions:##
#This defines a log-likelihood function for a spatiotemporal model, so the "optim" command can optimize it.
#The "par" argument is a vector of parameters. par[1] is the spatial lag coefficient. par[2] is the error variance coefficient. Remaining terms are regression coefficients.
#The "X" argument is a matrix of predictor values. Remember to include a constant vector to obtain an intercept. Users who want a temporal lag in the model need to include the term as a predictor in their "X" matrix.
#The "y" argument is a vector of dependent variable values. Be sure this vector is ordered in the same way as the rows of "X" are. Also be sure the length of y equals the number of rows of X.
#The "W" argument is a spatial weight matrix that defines the lags. This must be defined beforehand. See the file "laTimeSpace.R" for an example of how to create this.
#The "Time" argument refers to the number of waves of observations.
#The "N" argument refers to the number of cross-sectional units.

spatiotemporal <- function(par, X, y, W, Time, N) {
	n<-N*Time
	rho<-par[1]
	sigma2<-par[2]
	beta<-par[3:length(par)]
	ident<-diag(1,nrow=Time)
	nIdent<-diag(1,nrow=N)
	bigIdent<-ident%x%nIdent
	bigW<-(ident%x%W)
	first<-Time*(sum(log(1-rho*eigen(W)$values)))
    	loglikelihood <-first-(n/2)*log(2*pi)-(n/2)*log(sigma2)-((t(y-rho*bigW%*%y-X%*%beta)%*%(y-rho*bigW%*%y-X%*%beta))/(2*sigma2))
    	return(loglikelihood)
}
